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constructed from a pair of elementary projections and a single real parameter (5. 
For the standard application in optics, where the two projections implement Fourier 
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the "hybrid" form of Fienup's 1 input-output map for (3 = 1. Other values of (3 
are equally effective in retrieving phases but have no input-output counterparts. 
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1. Introduction 

In the standard phase retrieval problem one has available the modulus of the Fourier transform 
of an "object", usually in two or three dimensions, and a set of a priori constraints that a proper 
Fourier reconstruction of the object must satisfy. Usually one assumes that, up to translation and 
inversion of the object, the reconstruction is unique. In other words, there is essentially a unique 
set of phases that, when combined with the known Fourier modulus data, satisfy all the a priori 
constraints. It is in this sense that the unknown phases are said to be "retrieved" from the modulus 
data. 

Currently there are no practical algorithms for phase retrieval, that is, procedures where a 
solution is guaranteed at a computational cost that grows modestly with the size of the problem. 
Rather, a number of heuristic strategies have been developed that in many situations perform very 
well. A dominant theme of these practical methods is an iteration scheme which modifies the object 
in each cycle until both the Fourier moduli are correct and all the a priori constraints are satisfied. 
Although nothing has been proven about the rate of convergence of such schemes, the body of 
favorable empirical evidence in both imaging 2 and crystallographic 3 applications is substantial. 

Expressed in the most general terms, phase retrieval is the problem of finding an element in a 
large set that simultaneously has properties A, B, etc., given only the ability to find elements having 
these properties separately. If property A corresponds to the object having the correct Fourier 
modulus, property B might represent, for example, a support constraint. While it is very easy to 
construct objects having the correct Fourier modulus, or a prescribed support, there are no known 
methods for directly constructing an object (if one exists) that has both properties. Accepting this 
basic limitation, we can nevertheless attempt to build the desired construction using the simple 
constructions for properties A, B, etc. as elementary operations. 

A natural setting 4 for reconstructing an object is the ./V dimensional Euclidean vector space 
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E , whose components correspond to the values of N pixels, in an imaging application, or in the 
case of crystallography, a sampling of the electron density at N regular grid points in the unit 
cell. The main benefit of the Euclidean space is the ability to define a distance between objects. 
It is convenient to consider a complex Euclidean space, where unitary transformations such as the 
Fourier transform may act. Since distance is preserved by unitary transformations, distances can 
be measured in either the object or Fourier domains. 

With the ability to measure distance we can define more precisely the elementary operations 
for constructing objects with properties A, B, etc. 4 These are the "projections", which begin with 
some point p G E , and find a point tta(p) having property A, for example, and whose distance 
from p is a minimum. For most of the projections usually considered in imaging and crystallography, 
the point with minimum distance is essentially unique and can be computed in a time that grows 
no faster than iVlogiV. This is important, since we anticipate that many elementary projections 
will be required to arrive at the desired object having all the properties A, B, etc. 

The earliest application of projections to reconstruct an image from Fourier modulus and 
support constraints is the Gerchberg-Saxton 5 map: G = ir supp o 7r mo( j. When applied to an initial 
randomly chosen point, the idea is that successive iterates might simultaneously get closer, both to 
the subspace of E N having the correct Fourier modulus, as well as the subspace corresponding to 
the support constraint. Owing to the non-convexity 6 of the subspace of correct Fourier modulus, 
this map has the problem that it can encounter fixed points which have the correct support but fail 
to have the correct modulus. The Gerchberg-Saxton map in this situation moves the point from 
one constraint subspace to the other, and then back to the original point; a point common to both 
subspaces has not been identified. This phenomenon, potentially present for all projections where 
the corresponding subspace is non-convex, is known as "stagnation". 

To avoid stagnation, various generalizations 1 ' 6 of the Gerchberg-Saxton map have been pro- 
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posed. In imaging, the most successful of these have been schemes where the iterates are not confined 
to the constraint subspaces. In particular, objects realized as general linear combinations of other 
objects are generated that take advantage of the vector space nature of the search space. For these 
generalized maps it is usually the case that none of the iterates enroute to the solution satisfy any 
of the constraints. Fienup 1 has devised a family of maps motivated by ideas from control theory, 
of which the "hybrid input-output" map stands out as the most successful in a wide variety of 
imaging applications. The theoretical understanding of the hybrid input-output map, however, is 
very incomplete. There is not even a compelling argument that sets the "hybrid" version of the 
map apart from its less successful input-output relatives. Also, whereas the hybrid input-output 
map was developed for imaging applications with a support constraint, one would like to be able 
to exploit its excellent performance characteristics in applications such as crystallography, where 
support constraints are usually not available. 

This paper describes a general map formed by taking the difference of a pair of elementary 
projections. The advantage of such a "difference map" is that stagnation, in a strict sense, can be 
ruled out. By optimizing local convergence properties, the detailed form of the difference map is 
determined up to a single real parameter (3. Interestingly, for the special case of support constraints, 
the value (3=1 reproduces an input-output type of map, and it is exactly of the hybrid variety. 
Other values of /3, such as —1, give maps for which there are no input-output counterparts and 
yet whose performance in numerical experiments is also very good. Perhaps most significant is the 
versatility of the difference map with respect to the elementary projections from which it is built. 
In particular, all of the examples reconstructed in numerical experiments with support projection 
can, when implemented by the difference map, be reconstructed using an elementary projection 
onto a known distribution of object values. The set of object values, or histogram, is a form of a 
priori constraint that can be exploited in crystallography. The difference map thus provides a link 
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between the most successful phase retrieval method in imaging, the hybrid input-output map, and 
the crystallographic phase problem. 

2. Examples of elementary projections 

A projection, of arbitrary points in the N dimensional Euclidean space of objects, to points on a 
subspace, representing a particular constraint, is said to be elementary if it is (i) distance minimizing 
and unique, and (ii) easy to compute. In imaging and crystallography, "easy to compute" translates 
to the statement that the operation count scale essentially as N. Below is a collection of elementary 
projections that have been employed in various phase retrieval schemes: 

Fourier modulus. It is convenient to express all projections with respect to a common basis 
in E N , which we take to be the object domain. Because the Fourier modulus projection is more 
naturally expressed in the Fourier domain, we write 

TTmod = F~ ■ ^mod " F , (1) 

where T is the unitary transformation to the Fourier domain, and 7f mo d is the projection operator 
which acts componentwise, that is, on each pixel of the (discrete) Fourier transform. Geometrically, 
each complex component in the Fourier transform is mapped by 7r mo d to the nearest point on 
a circle having a prescribed radius (the Fourier modulus). Thus it follows that 7r moc j is distance 
minimizing and essentially unique (the exception being the measure-zero set of complex numbers 
with zero modulus). Since the computation of 7r mo d is linear in N, and the multiplication by T and 
requires N log ./V operations when using the FFT algorithm, Fourier modulus projection is easy 
to compute. Since circles are not convex, the subspace of objects satisfying the Fourier modulus 
constraints, C mo d, is also non-convex. 

Support. There is a unique, distance minimizing map onto the subspace C sup p of objects having 
a specified support S. If p n is the value of pixel n in the object domain, then support projection is 
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the map 

{p n if n e S 
( 2 ) 
if n i S. 

Since C supp is a linear subspace it is convex; projection computations require at most N operations. 
Positivity. For real-valued objects one can impose positivity. The unique, distance minimizing 
map to the corresponding subspace C pos , sets all negative pixels values to zero and leaves the 
positive pixel values unchanged. Like C supp , C pos is convex. Moreover, the projections 7r supp and 
^pos commute and can be combined into a single elementary projection. 

Histogram. Histogram projection 7 is complementary to support projection, if we regard the object 
as a function from the support S to a set of values, or "histogram", H. Here we define the object 
histogram H to be the set of N pixel values without regard to pixel position. For real- valued objects 
a distance minimizing map 7Thist ; onto the subspace of objects having histogram H, is unique and 
easily computed. One begins by sorting H = {hi, hi, . . .} and also finding an ordering of the pixels, 
n i — > o(n), such that their values, {p a m, |0 O (2)> • • •} are a ^ so sorted. Histogram projection is then 
given by the map 

TThist : Po{n) ^ h n . (3) 

It is straightforward to show that 7Thist is distance minimizing. Since N log N operations are required 
to sort N real numbers, histogram projection is easy to compute. For complex valued objects it 
is probably not possible to compute the projection to a specified complex histogram in order N 
operations, although approximate distance minimizing maps can be computed with this effort. 8 
The subspace of objects having a prescribed histogram, Chist) is the point set formed by applying 
all permutations to the point {h\, hi, . . .}, and as such is non-convex. 

Atomicity. The projection of greatest relevance to crystallography (and possibly astronomy) is the 
distance minimizing map to the set of objects consisting of a known number M of non-overlapping 
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atoms (or stars). It is not necessary to make the restriction to equal atoms, although this simplifies 
the computation of the projection. Our model of atoms will be objects of small support, usually 
just a 3 x 3 array of pixels (3 x 3 x 3 in three dimensions). The object values for each type of 
atom are also specified, and allowance is made for the possibility that the actual center of the atom 
may be located between the pixel centers (see Appendix for details). A unique, distance minimizing 
map is most easily constructed for the case of identical atoms with support on a single pixel. This 
situation may equivalently be characterized by its histogram: a set of M identical positive values 
p+, representing the atoms, and N — M zeros. Atom projection, in this case, corresponds simply to 
histogram projection: the M largest pixels in the object are set to p+ and the remainder set to zero. 
The difference between histogram and atom projection emerges when the support of each atom is 
larger than a single pixel; a graphical illustration is given in Figure 1. In general terms we can say 
that, whereas histogram projection completely neglects spatial correlations of the pixel values in 
the object domain, with atom projection one at least incorporates short range correlations. 

It is more difficult to construct atom projections for multi-pixel atoms, and consequently, to 
rigorously prove minimality. In practice it probably suffices to map the object reasonably near 
the true distance minimizing point on the subspace of atomic objects, C a tom- The algorithm used 
in the numerical experiments described in Section 7 begins by sorting the pixels by value and, 
beginning with the largest, associates them with the specified number of atoms while observing 
the restriction that each new atom's support does not overlap the support of previously identified 
atoms. The mapping is then fine-tuned by computing the centroid of the object values within each 
atom's support and setting the values on the support equal to that of an atom having that particular 
fractional location with respect to the pixels. More details are given in the Appendix. 

Apart from the minor complication due to non-overlapping supports, the subspace C atom is 
just the product of tori corresponding to the locations of atoms within the periodic crystal unit 
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(a) (b) 




Fig. 1. Examples of elementary projections, (a) Random positive density; (b) Fourier mod- 
ulus projection of (a); (c) histogram projection of (b); (d) atomicity projection of (b). 

cell (a non-convex set). This assumes that atom locations can move continuously between pixels. 
In practice this is only approximately true, since the support of each atom must move one pixel 
at a time. The discreteness of pixels also manifests itself in the relationship of C a tom to Chist, 
which should be one of containment: C at0 m C d^t- For this to be true we see that Chist must 
contain object values for all fractional displacements of an atom relative to the pixels. Clearly this 
is not true if we interpret Chist as the point set formed by taking all permutations of a given set 
of iV histogram values. The true histogram should include a continuum of values, corresponding 
to continuous translations of the atoms. For the practical implementation of histogram projection, 
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however, it suffices to sample just N representative values from the continuous distribution and use 
these as described above. 

Some observations of a general nature can be made if the constraint subspaces are smooth 
submanifolds of E . While this is the case for the Fourier modulus and support constraints, it 
applies approximately for atomicity since C a t m has a nearly smooth parametrization by continuous 
atomic positions on the d-dimensional torus. If one is seeking an object p 6 E N in the intersection 
of two smooth constraint spaces, say Ca and Cb, then just the dimensionalities of these spaces 
can provide useful information. For generic embeddings in E N , for example, one would not expect 
a unique solution if dimC^ + dim Cb > N. For phase retrieval (Ca = Cmod) to be well posed, 
this implies dimCs < N/2, that is, an object support that is smaller than half of the image, if 
Cb = Csupp; or a number of atoms M smaller than N/(2d), if Cb = C a tom- The converse of this is 
that when a generic embedding predicts an empty intersection, the fact that a solution is a priori 
known to exist implies it is unique. However, as the case of phase retrieval with support constraints 
in one dimension 9 shows, this form of reckoning can fail. Nevertheless, experiments (Sec. 7) with 
atomicity and histogram constraints suggest that nonuniqueness in the case of overdetermined 
constraints is the exception rather than the rule. 

3. The difference map 

Given two arbitrary projections m and tt 2 , we consider the "difference map" D: E N — > E N defined 



by 



D = 1 + PA 



(4) 



where (3 is a nonzero real parameter and 



A = 7Ti o f 2 - vr 2 o f 1 



(5) 



9 



is the difference of the two projection operators, each composed with a map fi'.E — ► E . The 
detailed form of the maps /j is secondary to the global behavior of the difference map and is 
discussed in the next Section. A fixed point of D, p*, is characterized by A(p*) = 0, or 

(tti o f 2 )(p*) = Pin2 = (tt 2 o fi)(p*) , (6) 

where pm2 is required to lie in the intersection of the corresponding constraint subspaces. If TT2 = 
7r m od, say, and tt\ represents some object domain constraint, then p\r\2 is a solution of that particular 
instance of the phase problem. We note that in general pm2 7^ P* ■ The standard approach for finding 
solutions is to begin with a randomly chosen point p(0) S E N , iterate the map D until a fixed point 
p* is found, and then use (6) to find pm2- Even if the solution pm2 is unique (up to translation 
and inversion), the set of fixed points is in general the large space given by 

(-7T1 ° .^r^Pire) n (7T 2 o /i) _1 (pm2) • (7) 

One practical consequence of this is that the fixed point object p* will appear to be contaminated 
with noise whose origin, ultimately, is the randomness of the starting point p(0). This is illustrated 
by the numerical experiments in Section 7. 

The progress of the iterates p(i) can be monitored by keeping a record of the norm of the 
differences 

ei = \\A(p(i))\\ , (8) 

where || • || denotes the Euclidean norm. The "error" has the geometrical interpretation as the 
currently achieved distance between the two constraint subspaces, C\ and C%. When this distance 
becomes sufficiently small, an object has been found that satisfies both sets of constraints. The 
behavior of ej is, in general, nonmonotonic. 

It is convenient to scale the constraint subspaces so that they lie on the sphere of unit norm. 
The magnitude of the error estimate is then related to the angular separation of the constraint 
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subspaces by 



e = 2 sin 



( 



0in2 
2 



) 



(9) 



While small values of 6m2 imply a solution is nearby, small values encountered during the early 
iterates are a cause for concern. This is because it is highly improbable that the random starting 
point p(0) is within a small angle of the solution, or more precisely, the set of fixed points (7). Small 
initial errors (angles) are more likely to be an indication that the object constraints are too weak 
for reliable phase retrieval, as for example an object support that occupies nearly half the image. 

4. Local convergence and traps 

The maps fi in the definition of the difference map (5) must be chosen with care if the set of fixed 
points is to be attractive. For example, choosing the identity map for both fi and fi does not 
give attractive fixed points. In order to limit the possibilities for these maps we use a geometric 
approach that relies on the elementary projections tt\ and 112- Let p G E N be the current iterate 
(object). Using projections we can obtain two additional points, iri(p) and ^(p). Pairs of points 
determine lines; in particular, 



is a general point, parametrized by a real number 73, on the line defined by p and vrj(p). Using (10) as 
a definition of the maps fi might be viewed as taking the first step beyond simply using identity maps 
for the fi. Indeed, the parameter choice 7; = — 1 reduces exactly to this case. However, we will see 
that we need to exploit the freedom associated with 7, to optimize the local convergence properties 
of the difference map. We also see that only the compositions 7Tj with fj, for i ^ j, make sense 
when considered locally. This is because locally the constraint subspace Cj may be approximated 
as an affine linear space, where vrjo/j = 7Tj, and hence nothing is gained by introducing the maps fi. 
Similarly, it is easily shown that nothing beyond the form (10) is gained by considering arbitrary 




(10) 
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points on the plane determined by all three points (p and its two projections, 7Ti(p) and ^(p))- 

The optimum values for 71 and 72 are determined by considering one iteration of the difference 
map for points p in the vicinity of a solution. Putting aside the case of histogram constraints, 
we assume that the constraint subspaces C\ and C2 are smooth, or nearly so. At a solution, C\ 
and Ci intersect and can be approximated as linear affine spaces. We will treat the case where 
the corresponding linear spaces X\ = C\ — C\ and X2 = C2 — C2 are orthogonal. Although this 
assumption may fail for certain pairs of projections, it prevails in the probabilistic sense of randomly 
oriented spaces satisfying dimXi + dim X2 < N in the limit N — > 00. If Y is the complement in 
E N of the span of X\ and X2, then a general point can be uniquely expressed as p = xi + X2 + y, 
where %i £ Xi and y €Y. The constraint spaces are now explicitly approximated as 

Ci = Xi + a 2 + h 

(11) 

C 2 = ai + X 2 + b 2 , 

where £ Xi and bi £ Y. Orthogonality of the spaces X\, X2 and Y implies the shortest element 
in C\ — C2 is 61 — 62 £ Y . Thus a solution (intersection) corresponds to 61 = 62, while 61 7^ 62 
represents a "trap", that is, a source of stagnation. The latter is illustrated by the action of the 
Gerchberg-Saxton map G = tx\o'K2- First note the formulas for projections of a general point (these 
make use of the distance minimizing property of ix\ and 7^): 

7ri (x\ + x 2 + y) = xi + a 2 + 61 

(12) 

7T2(xi + x 2 + y) =ai+x 2 + 6 2 . 
One application of G brings us to the fixed point 

G(x\ + x 2 + y) = ai + a 2 + 61 (13) 

and stagnation occurs. Subsequent applications of the elementary projections simply hop between 
a i + a 2 + 61 and 01 + 02 + 62, the two points on C\ and C2 with minimum separation. 
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The behavior of the difference map is quite different, as we now show. A straightforward 
calculation gives the result 

D( Xl +x 2 + y) = a l + a 2 + y + {l- P^fa - a x ) + (1 + fiyt)(x 2 - a 2 ) + f3(h - b 2 ) . (14) 

First consider the case b\ = b 2 = b, corresponding to a true intersection of the subspaces at the 
solution pi n2 = ai + a 2 + b. From (14) we see that subsequent iterates approach the fixed point 
p* = a\ + a 2 + y provided < (3^ 2 < 2 and —2 < /?7i < 0. This excludes 71 = 72 = —1, or identity 
maps, for the fa. Optimal convergence (one iteration) is achieved by the parameter values 72 = 
and 71 = — (3 ■ This choice also makes D maximally contractive; the rank of the Jacobian is then 
dimy. As explained in the previous Section, a fixed point of D guarantees a true solution and it is 
given by p m2 = (717 o f 2 ){p*) = a x + a 2 + b. 

Next we examine the situation for b± ^ b 2 , when the two constraint subspaces locally form a 
trap. With the parameters 7, set to their optimal values, subsequent iterates have the form 

D n {x l + x 2 + y)=ai + a 2 + y + nf3{b l -b 2 ) . (15) 

Rather than hopping between the two constraint subspaces we see that the iterates move away 
uniformly along the nearest separation axis, b\ — b 2 , with a rate determined by (3. In a sense, the 
map has recognized that there is no solution at hand locally and seeks one elsewhere. 

The local geometry of solutions and traps is rendered schematically in Figure 2. It should 
be understood that this view, and the analysis in this Section, is valid only when there exists a 
minimum local separation of the two constraint subspaces, <5 m ; n = — 62II) and (5 m ; n is small on 
the scale where these subspaces deviate from linearity. We summarize the results of the last two 
Sections with the final form of the difference map: 

D-.p^p + p^ [(l + zrVaOO-zrV] -7r 2 [(i-/rVi(p) + /3~V]) ■ (is) 
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(a) 



(b) 





Fig. 2. Constraint subspaccs (perpendicular rods) in the neighborhood of a solution (a) and 
a trap (b). The two circular disks in (a) represent points which project to the intersection 
of the constraint subspaces, pica-, under action of tt\ o / 2 and TT2 ° fi respectively; their 
intersection is the set of fixed points of the difference map. When the constraint subspaces 
do not intersect, as in (b), the action of the difference map is to move the iterates along the 
axis of minimum separation. 

We note that all nonzero values of (3 offer interesting choices, and that reversing the sign of (3 has 
the effect of interchanging the two projections. 

5. Relationship to Fienup's input-output maps 

From (16) we see that four elementary projections are computed for every iteration of the difference 
map. This number is reduced to two when (3 = ±1. In the case of support constraints, that is, 
7Ti = 7r sup p and iT2 = 7r m od> and the choice (3 = 1, one obtains the simple result: 



D 



f3=l ~ 



-^hybrid • P\ 




TTmod(p) 



Pn - vr mo d(/o)n if n S 



n 



if n £ S 



(17) 
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This is Fienup's hybrid input-output map for an object with support S and Fienup's parameter 
= 1. In his discussion of input-output maps, Fienup also considered the maps 



Only the hybrid map (17) can be obtained as a special case of the difference map. It is interesting 
that of the three input-output maps tested numerically by Fienup, the hybrid map outperformed the 
others and the optimum performance occurred for the parameter value ftp ~ 1. In these studies 1 
a fixed number of iterations of an input-output map was followed by several iterations of the 
Gerchberg-Saxton map in order to arrive at the solution. In view of the discussion in Section 3, 
concerning the distinction between the solution and the fixed point of the map, the Gerchberg- 
Saxton iterations can be avoided for (3y = 1. In this case the solution pm2 is obtained from the 
fixed point p* by an elementary projection (in effect a single Gerchberg-Saxton iteration): 



No such statement can be made for other values of /3f, for which there is no corresponding difference 
map. 

The choice /3 = — 1 in the difference map with support constraint is also relatively simple, but 
does not appear to have been studied previously. This has the effect of interchanging the roles of 
support and Fourier modulus projection. To express the resulting map most compactly we introduce 
a sign flipping operation on an object with support S: 




(18) 




(19) 



Pin2 = (tt 2 ° fi){p*) 



7Tmod(P*) 



(20) 




(21) 
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The P = — 1 counterpart of the hybrid input-output map is then given by 




(22) 



In this case the solution is obtained from the fixed point most directly by 



Pin2 = (tti ° h){p*) 



= 7T. 



supp 



(p1 



(23) 



Numerical experiments (Sec. 7) show that the difference map with (3 = — 1 is often just as effective 
in finding solutions as the choice (5 = 1. 

6. Atomicity projection and Sayre's equation 

Phase retrieval in crystallography has traditionally been performed purely in the Fourier domain, 
whereas the a priori constraint of greatest relevance, atomicity, is most directly expressed in the 
object domain. A key development in the history of crystallographic phase retrieval was Sayre's 
equation, 10 which permitted a simple translation from the object to the Fourier domain for the 
case of identical, non-overlapping atoms. In the object domain Sayre's equation reads 



where p x p denotes the nonlinear operation on the vector p 6 E which squares each component, 
and * is the discrete convolution operator. The object g is the point spread function corresponding 
to the finite atomic size a. To simplify the analysis we consider the continuum limit, where a is 
large on the scale of the unit spacing of pixels. Using r to denote pixel positions relative to some 
origin, then for a Gaussian atom in d dimensions, the choice 



p = g*{px p) 



(24) 




(25) 



is consistent with Sayre's equation when p = p^ , where 




(26) 
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We note that Sayre's equation "selects" not just the size of atoms, a, but also their normalization. 
The choice given above corresponds to 

IIP|| 2 = ^X>| 2 = 1 , (27) 

where M is the number of atoms. Sayre's equation is easily translated into the Fourier domain by 
interchanging componentwise multiplication and convolution. 

The earliest implementation of Sayre's equation in the Fourier domain was the "tangent for- 
mula," a mapping of phases given by 

T: arg p arg [g x (p * p)\ . (28) 

Since the Fourier modulus of p is fixed, T is a mapping of the subspace C mo d onto itself. Solutions 
to the phase problem for identical atoms must be fixed points of T, but the converse need not be 
true since equality of the phases does not preclude inequality of the moduli: 

\p\¥=\9*(p*p)\ ■ (29) 

In fact, a naive application of T is unstable to the trivial "uranium atom" solution, where all the 
phase angles are zero (or a translation thereof). 

A more sophisticated treatment 11 utilizes an objective function constructed from Sayre's equa- 
tion: 

V = \\\p-g*(pxp)\\ 2 ■ (30) 

Given some small e > which quantifies the error in the characterization of atomicity by Sayre's 
equation for the problem at hand, the region of E N satisfying V < e can be identified with the 
subspace of atomicity constraints, C at0 m- With the gradient of V, 

VV: p^ p-g*(px p)-2(g*p)x p + 2(g*g*(px p))x p , (31) 
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where g r = g- r , we can hope to find C a tom using the method of steepest descent. Within the usual 
framework of phase retrieval as practiced in crystallography, one would consider the map 

Smod = TTmod (1 ~ «W) . (32) 

In the limit a — ► + , the iterates of Smod fi° w i n the direction of decreasing V on the subspace of 
modulus constraints, C mo( j. Apart for being impractical, the evolution for a — > + is plagued by the 
problem of proliferating local minima (stagnation). Numerical experiments suggest that stagnation 
occurs over a range < a < a c , where at the upper limit only the most attractive minima, including 
the desired global one, are active. Thus by carefully tuning a near a c , the map *S* mo d can provide a 
reasonably practical phase retrieval scheme. A modification of the map for a = 1, which resembles 
the tangent formula when expressed in the Fourier domain, 



'mod| Q= i :argp -» arg 



g x (p*p) + 2(g x p) *p- 2 (\g\ 2 x (p*/5)) *p , (33) 



was used successfully in small molecule structure determination from x-ray data. 11 

Sayre's equation can also be used in conjunction with the difference map, where, by treating 
Fourier modulus and atomicity constraints separately, stagnation is avoided. Rather than S mo d, 
one considers 

Snorm = 7T norm O (1 - aW) , (34) 

where vr norm is the projection on the sphere with normalization (27). The fixed points of S^orm are 
atomic objects, without regard to Fourier modulus, and thus can be identified with C a t m- A good 
approximation to the projection operator is given by 

TTatom ^ ^Sayre = 'S'norm (^5) 

for large k. The local distance minimizing properties of vrs ayre , as required for convergence of the 
difference map, can also be checked. Formally this follows from the smooth form of VV, and the 
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structure of its eigenvalues when linearized about a V = fixed point. Since V = is the global 
minimum, all such fixed points are attractive or equivalently, the linearization of W has no negative 
eigenvalues. Zero eigenvalues will occur and correspond to continuous motion of the atoms, that 
is, motion along the subspace C a tom- For sufficiently small a, the linearization of Snorm is thus 
collapsing, on the eigenspaces with positive eigenvalues, and the identity on the local tangent space 
of Catom- By iterating, 5 n0 rm approaches a canonical (distance minimizing) projection operator. 

A handle on the parameters a and k is obtained from the stability analysis of a single atom 
having only a normalization degree of freedom: p = Xp^\ Moreover, when the total number of 
atoms M is large, we can neglect the action of 7r norm and obtain S'norm 

(A/ 1 )) = f a (X)p ( - 1 \ where 

f a : A A - a \x - A 2 + c(A 3 - A 2 )] , (36) 

and c = 2(4/3) d / 2 . The map f a has attractive fixed points A = 0, l,±oo, where the latter is the 
"uranium" instability. There is also a repulsive fixed point A* = c~ 1 which separates the "non- 
atom" (0) and "atom" (1) fixed points. For small a the flow from A* to either or 1 is slow, and 
many iterations k are needed. A reasonable criterion for selecting k is that the slope of /* at A* is 
greater than one by a significant factor, say 2: 

log2<log[/4(A*)] fc = A;log[l + a(l-c- 1 )] « ka{l - c~ l ) . (37) 

Figure 3 shows a plot of with a = 0.37, chosen so that (37) is an equality for d = 2 (c = 8/3). 
The uranium instability sets in shortly beyond the upper boundary of the A = 1 fixed point's basin 
of attraction, given by the largest root of 

= U [ Xo) 7 1 = 1 - a(cA 2 - A ) . (38) 
Ao — 1 

If we view Ao as the largest starting value that does not encounter the uranium instability, then 

a< 7XTY d ■ (39 > 
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Fig. 3. Plot of the map f a , Eq. (36), composed with itself three times for a = 0.37 and 
c = 8/3. Points of intersection with the straight line give the three fixed points A = 0, c _1 , 1; 
the divergent behavior near A = — 1 and A = 2 corresponds to the "uranium" instability. 

When the Sayre projection (35) is applied to an object p comprising a large number of atoms M, 
the starting normalizations Ao of the different atoms are drawn from a normal distribution. If the 
largest of these normalizations is in conflict with (39) the projection will fail. The maximum of M 
elements drawn from a normal distribution has the extreme value distribution 12 with mean that 
grows as ^log M. Inequality (39) thus implies that a be decreased roughly as the logarithm of the 
number of atoms. 

With appropriately chosen k and a, the difference map implementation of Sayre's equation 
does not suffer from stagnation, as is common with tangent formula based schemes. Comparisons 
with the more direct approach to atomicity projection, as described in Section 2, are given in the 
next Section. 
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7. Numerical experiments 

Apart from the case of support constraints in one dimension, 9 the role of dimensionality in phase 
retrieval seems to be minor. This is the conclusion of experiments performed in one, two and 
three dimensions using the difference map with histogram and atomicity constraints. The length of 
computations (difference map iterations) required to locate M atoms in the crystal unit cell, or M 
stars in an image of the sky, or M spikes in a time series, are empirically very comparable, given equal 
numbers of Fourier moduli in the appropriate dimension. Other properties of the object, for example 
compactness, appear to be far more important in determining the computational complexity of 
phase retrieval. Given this ambivalence, most of the experiments below are two dimensional in 
conformity with the dimensionality of the print medium. 

A uniform set of normalization and initialization conventions were used in all the experiments. 
All objects p (points in E N ) were normalized so that ||p|| = 1. When the Fourier modulus and 
histogram is normalized with the same convention, normalization is not required during the course 
of the iterations. With support and positivity constraints, normalization was applied after every 
projection. The components of the initial object p(0) were generated by a pseudo-random num- 
ber generator and then normalized. All the experiments used Fourier modulus projection as the 
projection 1T2. 

Figure 4 shows the 192 x 192 pixel image and corresponding Fourier modulus used in the 
first experiments. The difference map is first demonstrated using support and positivity constraints 
(-7Ti = 7T SU p P o 7r pos ) and the value (3 = 1, for which it reproduces Fienup's hybrid input-output map. 
A circle, measuring 112 pixels in diameter, was used as support. After only about 10 iterations from 
the random start an object comprising a circular ring of dark contrast was formed. Subsequently, as 
seen in the series of plateaus in the error plot shown in Figure 5(b), the various characteristics of the 
object were refined. The 14-fold modulation of the outer ring appears after about 1000 iterations; 
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Fig. 4. Object (a) and corresponding Fourier modulus (b) used in subsequent phase retrieval 
experiments. 
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Fig. 5. Phase retrieval using the difference map with support and positivity constraints, 
(a) Final iterate (fixed point of the map); (b) behavior of the error e$ with iteration i. 
The random gray contrast in (a) is removed, and the original object (Fig. 4(a)) is almost 
perfectly recovered, by a single application of Fourier modulus projection. 



22 




Fig. 6. Same as Figure 5 but with histogram constraints. Far fewer iterations are required 
to recover the same object. 

the text is revealed after another 3000 iterations, etc. After 5000 iterations the magnitude of the 
difference between successive iterates (as given quantitatively in the error plot) is so small that 
the object has arrived at the essentially static fixed point shown in Figure 5(a). The best estimate 
of the solution, pm2) is obtained by a single application of 7r moc j to the last iterate and is almost 
indistinguishable from the true object (Fig. 4(a)). We note that the random gray contrast in the 
fixed point image (Fig. 5(a)) is a manifestation of the fact that fixed points of the difference map 
live in a large space, the particular point chosen being subject to the randomness in p(0). Very 
similar results (not shown) are obtained with (3 = —1, whereas performance degrades both for 
small and large (5 of either sign. 

Histogram constraints are much more effective in retrieving the phases of the object in Figure 
4(a). With 7Ti = 7Thist and j3 = 1, the object is perfectly recovered in only about 100 iterations, 
as shown in Figure 6. Again, the mottled background in the fixed point object, Figure 6(a), is 
completely eliminated by a single Fourier modulus projection. Results with (5 = —1 are again very 
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(a) (b) 




Fig. 7. Plots of sorted pixel values, or "histograms." (a) Cornell seal, Fig. 4(a);(b) random 
disk, Fig. 8(a). 

similar, although there are qualitative differences in the appearance of the iterates. The larger 
initial error using histogram projection also indicates that the corresponding constraint subspace is 
significantly smaller than the subspace of support constraints used in the previous experiment. We 
conclude that for this particular object the histogram appears to contain more information than 
knowledge of the object support. 

The "histogram" (sorted list of 192 x 192 object values) used to implement histogram projection 
for the previous example is shown plotted in Figure 7(a). Because of the quirks in the histogram, 
peculiar to this particular object, the histogram data cannot be viewed as a valid form of a priori 
knowledge. An example that comes closer to a valid application of histogram constraints is the 
object shown in Figure 8(a). This object has the same circular support as the previous object. 
However, its pixel values were assigned random numbers from a uniform distribution; the resulting 
histogram is shown in Figure 7(b). Although the progress toward the solution was slower initially, 
the detailed pattern of random pixel values of this object was perfectly recovered after about 600 
iterations (Fig. 8(b)). 

Histogram constraints apply ideally in situations where the object is composed of individual 
elements, e.g. atoms, whose values are known in detail even if the location of each element is not. If 
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Fig. 8. Phase retrieval applied to a random disk (a) using histogram constraints; (b) plot of 
the error. The recovered object (not shown) not only had the correct circular outline but 
reproduced in detail the pixel values within the disk. 

we refer to objects with this property generically as "atomic" , then the histogram of atomic objects 
is specified uniquely by the number and distribution of their constituent "atoms" . For simplicity, all 
the experiments below were performed with M identical atoms, each atom being an approximation 
to a Gaussian on a small number of pixels. In the two dimensional experiments the support of 
each atom was a 3 x 3 array of pixels; the continuously variable center of the atom with respect 
to the central pixel then gives rise to a continuous distribution of values for the histogram. Details 
regarding the construction of approximate Gaussians for a general choice of atomic support are 
given in the Appendix. 

Objects having a specified number of atoms M were generated by applying atomicity projection 
to a random image, 



To simulate the effects of clustering, additional power was optionally applied to the low spatial 




(40) 
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frequency components of p ra ,nd- Specifically, the Fourier modulus of /9 ran d was specified using the 
formula 

IPrandl « = 1 m ; o ' (^1) 

q m +% 

where \q\ is the magnitude of the Fourier wavevector and qo = 27r/£ is the characteristic wavevector 
corresponding to a clustering length scale £. The limit £ — ► reproduces a white power spectrum 
and uncorrelated atomic positions in p a tom (apart from avoided overlaps). For finite £ the atoms 
in Patom are clustered in groups with linear dimension of order £. The phases of prand were drawn 
at random from a uniform distribution, giving the maximally random atomic object possible. We 
note that the Fourier modulus after atomicity projection, |/0 atO m|> is no longer a smooth function of 
q although its behavior at small \q\ resembles Eq. (41) when locally averaged over wavevector. At 
large \q\ the rapid decay of |/5 a tom| with \q\ reflects the approximate Gaussian Fourier transform of 
each individual atom. 

Figure 9(a) shows an example of 60 uncorrelated atoms (£ — > 0) on 128 x 128 pixels. The 
adjoining error plots document the successful phase retrieval for this object by difference map 
iterations with three different choices for n\: (b) histogram projection and (3 = 1, (c) atomicity 
projection with 3x3 pixel Gaussian atoms and (3 = 0.5, (d) atomicity projection using Sayre's 
equation (tti = 7rs ayr e> Eq. (35), k = 3, a = 0.37) and (3 = 1. The error plots suggest that 
the phases are retrieved by a more nearly random process than is typically the case for compact 
objects, such as the ones considered earlier. A large number of experiments support the view that 
with uncorrelated atoms the search for the solution is a deterministic but nevertheless random 
walk through a relatively homogeneous "landscape" . The step size in this walk, as measured by the 
angle 6in2, is very nearly constant and only shrinks significantly when, apparently by accident, the 
walk arrives at the attractive basin of the solution. That there is no manifest progress toward the 
solution until the very end of the process is confirmed by snapshots of the object along the way. 
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Fig. 9. (a) Object comprising 60 equal "atoms"; error plots for three different implementa- 
tions of a priori information: (b) histogram constraint, (c) 3 x 3 pixel atoms, (d) atomicity 
implemented by Sayre's equation. 
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Fig. 10. Comparison of phase retrieval successes after 100 iterations of the difference map for 
histogram projection, 3x3 pixel atomicity projection, and atomicity projection via Sayre's 
equation. The horizontal axis is the difference map parameter f3. See text for more details. 

In contrast, the reconstruction using histogram information of the compact objects in the previous 
experiments proceeded through definite stages of refinement. Although the Sayre projection, Fig. 
9(d), found the solution relatively quickly, experiments with larger numbers of atoms displayed 
the same random walk behavior exhibited by the other projections. We also note that because the 
Sayre projection is much more demanding computationally than atomicity projection in the object 
domain (7r a tom), the performance of the two forms of atomicity projection is not as different as the 
error plots imply. Object domain atomicity projection has the added advantage that multiple kinds 
of atoms are easily accommodated. 

The random walk behavior of the difference map for uncorrelated atoms is also observed with 
other values of (3. A quantitative comparison that is free of object idiosyncrasies is achieved by 
averaging results over the ensemble of atoms defined above. The results of a study with uncorrelated 
(£ — ► 0) atoms on 64 x 64 pixel images is shown in Figure 10. For each value of (3 shown, and 
three choices for tt\ (7Thi s t, vr at0 m and vrsayrc)) 20 experiments were performed, each with a different 
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Fig. 11. Phase retrieval becomes easier when the atoms are clustered, as in this example (a) 
comprising 200 atoms; error plot for the case of 3 x 3 pixel atomicity projection. 

realization of atoms and 100 iterations of the difference map. In the studies of 7Thi s t and 7r a t m 
systems of 30 atoms were used and the success rate has peaks at both positive and negative (3. 
Since for this number of atoms 7rs a yre gave a 100% success rate over a large range of /?, the number 
of atoms was increased to 40. The success rate for this projection has a broader peak at larger (3 
and virtually zero success at negative (3 due to the "uranium" instability. We note that 7rs ayre has 
the additional parameters k and a that can be tuned to improve performance. This study used 
k = 3 and a = 0.37. 

The number of difference map iterations required for phase retrieval grows rapidly with the 
number of atoms and is the subject of a future study. Here we make the observation that any 
degree of clustering in the atomic positions has a profound effect, apparently reducing the number 
of iterations by orders of magnitude. Figure 11(a) shows an example with 200 atoms generated 
with a clustering length scale £ = 50 on 128 x 128 pixels. Without clustering, a system of this many 
atoms could not be solved using the difference map in a reasonable time. The effect of the clustering 
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Fig. 12. Phase retrieval for an object comprising 60 atoms in three dimensions, (a) The first 
16 layers, arranged lexicographically, of the 32 x 32 x 32 voxel array; (b) error plot for the 
difference map with histogram projection. 

is to greatly reduce the search space, in a manner not unlike what was observed for the compact 
objects (Figs. 4(a) and 8(a)) studied earlier, where enhanced power in the low spatial frequencies 
helps locate the "support" of the clusters. Snapshots of the iterates show a relatively slow variation 
of low frequency features accompanied by rapid fluctuations in the atom positions. The example 
shown was solved in about 600 iterations using ir\ = 7r atom and (3 = 0.5. The fluctuations in the 
error (Fig. 11(b)), after the solution has been found, correspond to random translations of the 
entire object. 

For completeness we include two experiments with atomic objects in dimensions other than two. 
As already mentioned, the performance of the difference map was found to be indistinguishable from 
studies in two dimensions having the same total number of pixels and atoms, all other attributes 
(clustering, atomic size, etc.) being the same. Figure 12 shows an experiment with 60 uncorrelated 
3x3x3 voxel atoms on a 32 x 32 x 32 grid; 16 layers of the object are arranged in Fig. 12(a). The 
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difference map using histogram projection and (3=1 found the solution in about 250 iterations. An 
experiment in one dimension is documented in Figure 13. The plot of the object values, Fig. 13(a), 
shows 60 uncorrelated "atoms" , each with support on 3 adjacent grid points; the size of the object is 
2 14 = 16384, the same as the two dimensional example in Figure 9(a). Figure 13(b) shows the fixed 
point of the difference map, again with histogram projection and (3 = 1. One application of Fourier 
modulus projection to the object in Fig. 13(b) perfectly reproduces a translated and inverted copy 
of the original object, Fig. 13(a). We see that the number of iterations required to find the solution 
(Fig. 13(c)), about 2000, is very similar to the result obtained in the two dimensional experiment, 
Fig. 9(b). 

8. Conclusions 

We are still far from being able to claim that an algorithm for phase retrieval is at hand. Like other 
iterative schemes currently in use, the difference map considered here falls short of a true algorithm 
in that there are no useful bounds on the number of iterations required to find the solution (a fixed 
point). The analysis given in Section 3 only guarantees that every fixed point corresponds to a 
solution (by Eq. (6)). Fixed points of dynamical systems are just the simplest form of an attractor; 
deterministic systems are also known to have "strange attractors", complex spaces upon which 
iterates of the map can become trapped. From numerical experiments there is some evidence that 
the existence of strange attractors, in competition with the simple fixed points, is the main failure 
mechanism (stagnation) of the difference map when (3 is too small. For small j3 the iterates exhibit 
chaotic fluctuations while preserving certain gross features in the object. Strange attractors can 
be tolerated when their number is small: one need only restart the iterations, hopefully within the 
basin of attraction of a simple fixed point. This is not an option when strange attractors proliferate, 
apparently for small {3, and possibly also in situations where the phase retrieval problem is close to 
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Fig. 13. Phase retrieval for an object comprising 60 "atoms" in one dimension, (a) Plot of 
object values; (b) fixed-point of difference map with histogram projection; (c) error plot. 
A single application of Fourier modulus projection to (b) perfectly reproduces (a), up to 
inversion and translation (not shown). 
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being underdetermined. 

Competing strange attractors appear to be absent when (3 is not small and the phase retrieval 
problem is well posed (overdetermined) . In this regime the numerical evidence suggests that the 
solution will always be found given enough iterations. By construction (Sec. 4), the difference map 
is highly contractive and we are led to believe that the fraction of "phase space" actually explored 
by the map is extremely small. 

Some insights about the nature of the space being explored can be gained from the study of 
phase retrieval experiments where the data has been fabricated to remove all fixed points (solutions) 
while preserving the general characteristics of the map. One such fabrication, for example, is to 
form two different objects having the same histogram and attempt phase retrieval using their 
averaged Fourier moduli. It is highly unlikely that phases exist which can be combined with this 
fabricated modulus data and yield an object having the required histogram. On the other hand, 
because the statistical properties of the fabricated data is very similar to that of either of the two 
genuine objects, the qualitative behavior of the difference map should be unchanged. The only real 
difference between the genuine and fabricated data is that the difference map iterates for the former 
will eventually find a fixed point, while for the latter the iterates will continue exploring indefinitely 
a space that we expect to share all the essential characteristics of the space being explored in the 
case of the genuine objects. This space is itself a strange attractor and the complexity of iterative 
phase retrieval is perhaps best quantified by some measure of its size. 

A "phase portrait" of a strange attractor, obtained with fabricated data of the kind described 
above, is shown for three values of (3 in Figure 14. The two genuine objects were one dimensional 
with identical two-valued histograms: distinct sequences of 16 zeros and 16 ones. Iterates were 
obtained for the difference map with histogram projection. Two dimensional images of the attractor 
were generated by first projecting each iterate into the four-torus corresponding to the first four 
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Fig. 14. Evolution of the difference map attractor with decreasing /?: (a) 0.7, (b) 0.6, (c) 0.5. 



(unrestricted) phase angles of the Fourier transform. Selecting only those iterates whose second 
pair of phase angles were both close to zero, the first pair of angles were plotted as a point in the 
plane. The images are thus projections of codimension-2 sections of the attractor. Figure 14 shows 
the evolution of the attractor from a relatively uniform probability distribution for [3 = 0.7 to a 
highly nonuniform distribution for j3 = 0.5. The collapse of the search space with decreasing (3 is 
consistent with a corresponding improvement in phase retrieval performance (for non-fabricated 
data). The absence of inversion symmetry in the image for (3 = 0.5 indicates that the threshold to 
the regime of multiple attractors has been crossed. 

The study of iterative phase retrieval in the context of discrete dynamical systems may provide 
new insights for improving solution strategies and even a handle on a comprehensive theory of the 
complexity of algorithms. A meaningful measure of the size of the search space, for example, might 
be given in terms of the dimensionality of the associated attractor. It would be interesting to study 
not just the evolution of complexity (e.g. attractor dimensionality) with parameters defining the 
iterative map (/3), but also with respect to quantifiable attributes of the object (clustering). 

The application of phase retrieval which currently poses the greatest challenge is macromolecu- 
lar crystallography. Recently the Shake and Bake (SnB) method 14 has succeeded in retrieving phases 
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for small proteins without the benefit of the additional MIR or MAD data normally required for 
protein structure determination. SnB 13 is iterative and uses object domain atomicity projection 
and objective function minimization as its two elementary operations. On general grounds, how- 
ever, one can argue that the SnB method may not have realized its full potential. First, the two 
operations are implemented in the alternating Gerchberg-Saxton fashion which in most other appli- 
cations leads to stagnation. Second, the objective function used by SnB is based on an uncorrelated 
collection of identical point scatterers, not unlike the function V (Eq. (30)) derived from Sayre's 
equation. An objective function (or associated projection) that makes use of clustering, or the divi- 
sion of space into solvent and non-solvent regions, would add considerable a priori information and 
thereby reduce the size of the space being searched. An example of the improved performance with 
tighter object domain constraints was seen in the experiment above which compared histogram and 
atomicity projection (Figs. 9(b) and 9(c), respectively). 

9. Appendix 

A. Finitely sampled Gaussians 

An "atom" in d-dimensions with center ro is modelled as the Gaussian 



We are interested in approximating f ona small set of pixels S = {si, S2, ■ ■ •} on a <i-dimensional 
cubic grid with unit spacing. For a given choice of "atomic support" S, there is an optimum Gaussian 
width a that we aim to determine. 

The support S is defined in a translationally invariant way relative to the atom center ro, 
or rather, the grid point po closest to tq. For randomly placed atoms we expect the fractional 




(42) 



with normalization 




(43) 
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translation t = tq — po to be distributed uniformly in the cube T = (— ^, ^) d . 

Let {^(t) | s £ S} be the approximation on S of the Gaussian \P with fractional center t £ T. 
The normalization 

£l^(*)| 2 = l (44) 

applies for each t £ T. Also, for each t € T, the values ^(i) are determined by minimizing the 
Euclidean distance between ^ and the restriction of to S 1 : 



<5(t) = £|§ s (t)-*( S ,;po + 0| 2 • (45) 
Minimizing 5{t) subject to Eq. (44) gives 

*.M- , * {SM + t) , (46) 



/E.6S + ')l 2 

and the squared Euclidean distance 



*(*)= ^El*( s ^o + i)| 2 -lJ • (47) 

Finally, the optimal Gaussian width a is determined by minimizing the average of 5{t) over fractional 
translations: 

<Wc = / S(t)d d t . (48) 

The minimization of <5 ave with respect to a can be carried out numerically on a suitably fine grid of 
fractional translations t. Results for various choices of support in d = 1, 2 and 3 are given in Table 
1. A natural choice for S is the set of pixels within a distance R = 1, y/2, \/3, ... of the origin. We 
expect the width a to increase and the distance <5 av c to decrease as R increases. 

The order of choosing the atomic support and corresponding Gaussian width a is usually 
reversed in actual applications. For example, in crystallography the Gaussian decay of intensity 
with diffraction wavevector q (due to atomic size, disorder, etc.) fixes the value of a. For M identical 
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Table 1 . Widths a and errors 5 avc of finitely sampled Gaussians in dimension d for supports 
within a distance R of the origin. 



Gaussian atoms at positions r m we have 
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(49) 



where 



M 



(50) 



m=l 



has essentially a white spectrum for random r m . A good estimate of a is therefore given directly in 
terms of the expectation value of \q\ 2 with respect to the intensity data \p q \ 2 - 



a 



1EW 

D EM 2 



(51) 



After determining a one would then consult a table such as Table 1 to discover the atomic support 
S for which this a is optimal, or close to optimal. 
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B. Atomicity projection 

We assume that finitely sampled Gaussians {^(t) | s £ S} have been determined for the problem 
at hand. These are precomputed on the chosen support S with a suitably fine grid of fractional 
translations t. An object comprising M identical atoms will have support S+{pi,p2, ■ ■ ■ ,Pm}, where 
the p m are atomic support centers on the integer grid. Normally one imposes a non-overlapping 
condition such as 

Pm-Pn^S-S (m^n). (52) 

Atomicity projection of an arbitrary (real valued) object p is accomplished in three steps: 
Obtain support centers. To identify positions in p having a large overlap with a Gaussian atom, 
g, one forms the convolution p' = g * p. The pixel values of p' are then sorted and their locations 
are appended, beginning with the largest, to the list of atomic support centers p m . Each potential 
new center must satisfy a separation condition such as (52). 

Obtain fractional positions. Within each atomic support S + po the actual atomic position may 
have a fractional translation t as described above. Since the finitely sampled Gaussian {^f s (t) \ s S 
S + Po} is an approximation to the true Gaussian ^>(r, po +*)> its centroid (on S + po) will be close 
to the true Gaussian center po + t. To find t one therefore obtains the centroid of p' restricted to 
S+po. If p' is far from being atomic it may happen that the t determined via the centroid is outside 
the cube T of allowed fractional translations. It is then necessary to obtain the translation if on 
the boundary of T that is closest to t. 

Object synthesis The final step is to add up the individual atomic objects. If atomic supports 
are allowed to overlap, the result must be normalized. 

Although this projection algorithm is probably not distance minimizing for a general object 
p, it is nearly so in the situation which matters most, that is, when p is already nearly atomic. 
We recall that the distance minimizing property of projections was required only in the analysis of 
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local convergence (Sec. 4). 
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